Valence Bond Solids and Their Quantum Melting in Hard-Core Bosons on the 

Kagome Lattice 
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Using large scale quantum Monte Carlo simulations and dual vortex theory we analyze the ground 
state phase diagram of hard-core bosons on the kagome lattice with nearest neighbor repulsion. In 
contrast to the case of a triangular lattice, no supersolid emerges for strong interactions. While a 
uniform superfluid prevails at half-filling, two novel solid phases emerge at densities p = 1/3 and 
p = 2/3. These solids exhibit an only partial ordering of the bosonic density, allowing for local 
resonances on a subset of hexagons of the kagome lattice. We provide evidence for a weakly first- 
order phase transition at the quantum melting point between these solid phases and the superfluid. 

PACS numbers: 05.30.Jp, 75.10.Jm, 75.40.Mg, 71.27.-fa 



Current interest in microscopic models of frustrated 
quantum systems stems largely from the search for ex- 
otic quantum phases and spin liquid states. In general, 
geometric frustration tends to destabilize quasi-classical 
order, possibly allowing for nontrivial quantum states 
and novel critical phenomena to emerge in such systems. 
One intriguing approach addresses classically frustrated 
(Ising) models perturbed by quantum (off-diagonal) in- 
teractions Q. The behavior of classically disordered, 
degenerate ground state manifolds upon application of 
a U(l) symmetric perturbation (e.g. ferromagnetic ex- 
change) is of special interest as experimental advances in 
the construction and control of atomic gases in optical 
lattices have opened up the possibility of designing such 
Hamiltonians for ultracold bosons. In particular, it has 
recently been shown how an optical kagome lattice can 
be constructed using a triple laser beam design Q , which 
could permit access to parameter regions of interest in the 
search for exotic quantum phenomena. 

In this paper we consider a model of bosons on the 
kagome lattice in the strongly interacting regime, cor- 
responding to the hard-core limit of the Bose-Hubbard 
Hamiltonian discussed in Ref. 0], 
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where bj (fe;) creates (destroys) a particle on site i, t > 
denotes the nearest-neighbor hopping, > is the 
nearest-neighbor repulsion, and fi is the chemical poten- 
tial. This model can also be mapped onto the spin- 1/2 
XXZ model allowing for an interpretation of our re- 
sults in terms of both bosons and quantum spins. We re- 
port results on the ground state phase diagram obtained 
from a combined analysis of large scale quantum Monte 
Carlo (QMC) simulations using the stochastic series ex- 
pansion technique yUsIl and phenomenological dual vor- 
tex theory (DVT) jaQ- We find that, in contrast to 



FIG. 1: (color online). Ground state phase diagram of hard- 
core bosons on the kagome lattice (inset). The primitive vec- 
tors ai and a2 are constrained on a (periodic) torus spanned 
by Li = m X ai and L2 = 712 x a2 (where ai = a2 = 2). The 
circles illustrate the subset of hexagons with a resonating bo- 
son occupation of three bosons per hexagon in the p — 2/3 
solid. The remaining bosons localize to form a solid backbone 
on the sites that do not belong to any of these hexagons. 



previous theoretical expectations, a uniform superfluid 
persists at half- filling for all values of V/t. In addition, 
for fillings p — 1/3 and 2/3 we find evidence for valence 
bond solid (VBS) phases where bosons are delocalized 
around a subset of hexagons (see Fig.^l. We find that the 
quantum melting of both VBS phases into the superfluid 
occurs at weakly first-order quantum phase transitions. 

Past work on the ground state phase diagram of this 
model has been controversial and intriguing: Spin-wave 
calculations suggest that a supersolid state may emerge 
around half-filling (p = 1/2) at p = 2V for t/V < 0.5 H]. 
However, these results are not conclusive since strong 
quantum fluctuations may destroy the long-range order 
assumed within mean-field theory [3] ■ More recently, con- 
sideration of the large classical degeneracy at i = 
has led to the proposal of several exotic Mott-insulating 
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FIG. 2: Contour plots of the structure factor (left panel) and 
the static susceptibility (right panel) in the VBS close to the 
phase boundary for ii.IV = 11/12, tjV = 1/8, and T = i/24. 
The BZ of the underlying triangular lattice is hexagonal (see 
Fig- 01 with corners at Q = (27r/3,0) and symmetry related 
momenta. In both panels the axes range from — 27r to 27r, and 
ferromagnetic peaks at reciprocal lattice vectors have been 
suppressed. 



states (e.g. VBSs or disordered quantum liquids) at half- 
filling 10]. Very little work has been done to elucidate 
the nature of the phase diagram away from half-filling. 

Using QMC simulations, we have obtained the phase 
diagram of illustrated in Fig. ^ The lattice is 
empty for [i < —At and completely filled for > 
A{t -\- V). For large values of t/V , the bosons are su- 
perfluid, with a finite value of the superfluid density 
Ps, which we measure through winding number fluc- 
tuations Wai 2 ill each of the lattice directions as 

= ((W^aJ + {Wl^))l{2pt), where /3 is the inverse tem- 
perature. In agreement with mean-field theory, we find 
that two solid phases with rational fillings p — l/i and 
2/3 emerge at smaller t/V . Both are characterized by 
finite values (in the thermodynamic limit) of the den- 
sity structure factor per site, S{q)/N — (pqTPqr) and 
the static susceptibility per site, = (/ dTp^Tp\i:))^ 

where p^r = (1/iV) Pir exp(iqri) and pir is the boson 
density at site i and imaginary time r, at wave vectors 
q = Q = (27r/3, 0). This can be seen in Fig. 12 which 
shows peaks at the corners of the Brillouin zone (BZ) in 
the solid phase that are absent in the superfluid. 

Decreasing t/V at half-filling {p/V — 2), the value of 
Ps in the limit V ^ oo {t = \) takes on about 54% 
of its value at the XY point, Ps{V — Q) k, 0.55 (see 
Fig. 13) . This is in stark contrast to the triangular lat- 
tice case, where ps{V oo) approaches only about 4% 
of its value at the XY point [l3|- We do not observe 
evidence of long-range order in the density structure fac- 
tor, thus eliminating the possibility of a supersolid phase 
as found recently on the triangular lattice at small val- 
ues of t/V 0, 0, Q|. In addition, we do not observe 
any Bragg peaks in the bond-bond structure factor (de- 
fined below), precluding the existence of any VBS order. 
The persistence of the superfluid phase can be under- 
stood from a duality analysis of the boson problem in 
terms of vortices on the dual dice lattice, which interact 
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FIG. 3: (color online) . The superfluid density ps as a function 
of t/V at half-filling for p/V = 2 and T=t/5. The lattice 
linear dimension is denoted L (see Fig.0, so that the number 
of sites is = L X L X 3. Inset: The static structure factor 
S'(Q), superfiuid density ps, and boson density p as a function 
of p/V at t/V = 1/10, for L = 24 and T = t/24. 



with a dual magnetic flux of 2TTp/q = tt (for boson fill- 
ing p = p/q = 1/2) 0, 0- It has been shown that for 
the dice lattice at p = 1/2 the vortices undergo dynamic 
localization due to an Aharanov-Bohm caging effect 
that suppresses the condensation of vortices and leads to 
a persistence of superfiuidity 7]. 

At boson fillings p = 1/3 and 2/3, the model in the 
large V/t limit can be mapped to a quantum dimer model 
on the hexagonal lattice, where the occupied (unoccu- 
pied) sites at p = 1/3 {p = 2/3) correspond to the 
presence of dimers. In the classical limit {V/t = oo) 
all of the dimer coverings are degenerate. When V/t be- 
comes large but finite, the degeneracy is partially lifted, 
and the remaining dimer kinetic energy promotes the so- 
called plaquette order 0, 0| . In terms of bosons, for 
example at p = 2/3, every third site on the kagome lat- 
tice forms a solid backbone of occupied sites, whereas 
the remaining bosons locally resonate around every third 
hexagon, as shown in Fig. ^ DVT at p = 1/3 and 2/3 
predicts several possible valence bond order phases on 
the kagome lattice depending on the parameters of the 
effective field theory of the dual vortices. A set of mean- 
field states obtained by DVT are found to have the same 
spatial symmetry as the state described above. We ex- 
pect that the kinetic energy gain would prefer the phase 
with the resonating hexagons. Our QMC data on S{q) 
and x(q) is consistent with these expectations; in partic- 
ular, Bragg peaks are observed at the correct positions, 
as shown in Fig. 2. We find further QMC evidence for 
such resonances by measurements of the hexagon occupa- 
tion. Denoting by Pn the measured fraction of hexagons 
with n bosons, we obtain Po,i,2 < 10~^, P3 = 0.33(1), 
P4 0.38(2), P5 = 0.225(4), Pg = 0.052(4) at vari- 
ous points inside the p = 2/3 solid, which compare well 
to the expected values P^ of Pg^j = 0' = 5/12, 
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1/4, Pg = 1/12 for the resonating state. Fi- 
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FIG. 4: (color online). Finite-size scaling of the bond 
structure factor S'b(Q) in the VBS phase for pLjV = 11/12, 
tjV = 1/8, and T = t/12. Error bars are smaller than the 
symbol sizes, and the line shows a linear extrapolation to the 
thermodynamic limit. 
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FIG. 6: (color online). Data collapse of the superfluid density 
pB and static susceptibility x(Q) a-t p/F = 11/12 for /3/L — 
lit. 




FIG. 5: (color online). The correlation function Cb(ri — r^) 
between the bond indicated by large text (red online) and 
the other bonds for L = 24, njV = \\IV2, t/V = 1/8, and 
T = t/12. 



nally, we find sharp peaks in the bond structure factor 
and corresponding susceptibility at the corners of the 
kagome BZ. Fig. 0] shows the finite size scaling prop- 
erties of the structure factor S'b(q)/A^ — {B^t-B^t) 
q = Q, where B^r = (1/iV) Bq,^ exp(iqrQ,) summed 
over the bond index a connecting spins i and j, and 
Ba{i,j),T = ti^ibj + &i&])r is the ofF-diagonal bond oper- 
ator at imaginary time r. The non-zero intercept of the 
extrapolation in Fig. ^ exhibits long range order in the 
bond-bond correlations The real-space correlations, 
Chir^-rs) = ((^ / B^rdT-Ba){^ J BsrdT-Bo)), where 
Bq denotes the background bond strength, confirm the 
expected preponderance of resonating bonds on hexagons 
as illustrated in Fig. |S1 

As shown in Fig. [SJ "bow-tie" features appear in 5(q) 
and At low temperatures and close to qo = 

(0, 27r/\/3) and equivalent positions, these can be fit- 
ted to S'(qo + q) - {ql + K'^)/^/fT^, x{<lo + <l) - 
(q^ + K^)/{q'^ + K^), where is a component of q along 
the directions of bow-ties, and k is interpreted as an in- 
verse correlation length. The bow-ties are present in the 
solid phase and also in the superfluid phase near the tran- 
sition, albeit with larger values of k. Implications of this 
bow-tie structure will be discussed elsewhere (5j. 



In contrast to the case of a triangular lattice 'iS*!, we 
do not observe any obvious discontinuities that would in- 
dicate a first-order phase transition from the solid phases 
into the superfluid. An example is illustrated in Fig. 13 
for ps, S{Q) and p along a cut at fixed t/V from the su- 
perfluid into the p — 1/3 solid. Landau-Ginzburg- Wilson 
(LGW) theory forbids a generic continuous transition be- 
tween a solid and a superfluid, as different symmetries 
of Hb are broken in both phases. Possible explanations 
of our observation would thus include (i) a weakly first- 
order superfluid-solid transition, at which the correlation 
length stays finite but exceeds the linear size of the QMC 
cells; (ii) a narrow intermediate supersolid regime; or 
(iii) a non-LGW continuous transition, such as the re- 
cently proposed deconfined quantum criticality scenarios 
[20I l2l| . To discern between these possibilities we ob- 
tained detailed data, and performed a finite-size scaling 
(FSS) study and analysis of energy histograms over the 
melting transition region. 

We look first at the FSS behavior. For a two- 
dimensional system, the following FSS relations ap- 
ply for a continuous solid to superfluid quantum phase 
transition. The superfluid density scales as ps = 
L-^Fp^L^/^iKc - K),I3/L^), where L denotes the lin- 
ear system size, z the dynamical critical exponent, v the 
correlation length exponent, Kc — K the distance to the 
critical point in terms of the control parameter K , such 
as t/V or p,/V , and Fp^ the corresponding scaling func- 
tion. Similarly, S'(Q) ^ L^-''-'^ FsiL^'" {K^ - K), (3 / L^) 
and x(Q) = L'^-'^F^{L^/''{K,-K),P/U'), where is the 
anomalous exponent. From these scaling relations, given 
appropriate values for {t/V)c and the critical exponents, 
QMC data for different system sizes should follow uni- 
versal curves Fp_,{-,A) and F^{-.A), if the transition is 
continuous. 

QMC data was obtained over the melting transition, 
focusing on a fixed value of fi/V — 11/12, which is 
close to the largest extend of the p = 1/3 solid. As 
discussed below, the scaling of the QMC data appears 
consistent with z = 1, so that most simulations were 
performed for fixed aspect ratios A = P/L. As shown 
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FIG. 7: (color online). Distribution (arbitrary units) of the 
kinetic energy close to the transition point for various sys- 
tem sizes, [ilV = 11/12, tjV = 1/7.796, and T = f/60. As 
the system size increases, the double-peaked features become 
more pronounced. 

in Fig. El the data for p^L appears to collapse very well 
for {tlV)c = 0.12821(2), z = \, wA v = 0.43. The 
data for x(Q) also collapses well for (tlV)c = 0.12827(4), 
V = 0.45, and t] = —0.50(15). The two independent esti- 
mates for {t/V)c are sufficiently close in value to confirm 
a direct transition and make scenario (ii) seem unlikely. 
This is also consistent with the discovery of similar ly val- 
ues on both sides of the transition. Based on the value 
of 77 = —0.50(15) we also find the scaling of S'(Q) to be 
consistent with a dynamical exponent z — 1.0(2). 

Even though the FSS analysis appears consistent with 
a continuous quantum critical point in most respects, a 
negative value of 77 is unusual, and one may be inclined 
to interpret this as an indication of very weak first-order 
behavior. We explore this possibility in further detail by 
studying histograms of physical quantities in the transi- 
tion region. In Fig. |3 one can see double-peaked struc- 
tures developing in histograms of the kinetic energy over 
the melting transition for sufficiently large system sizes. 
We believe that this additional observation provides co- 
gent evidence for scenario (i), i.e. the melting transition 
studied here is indeed very weakly first-order. 

In conclusion, we have studied hard-core bosons on the 
kagome lattice with nearest-neighbor repulsion. At half- 
filling the model remains in a uniform superfluid phase 
for all values of V/t, and no supersolid emerges. For the 
solid phases at filling fractions 1/3 and 2/3, we find ev- 
idence for an exotic insulator with partially delocalized 
bosons occurring on a six-site hexagonal structure on the 
kagome lattice. Although the superfluid-solid melting 
transition appears naively to have the scaling properties 
of an unusual continuous quantum critical point, we find 
clear indication of weak first-order behavior from double- 
peaked histograms of the kinetic energy. The apparent 
weakness of this transition is in stark contrast to exam- 
ples of strong first-order quantum melting transitions in 
related models ■ An understanding of this contrasting 
behavior may have critical importance in the search for 
unconventional quantum criticality in this class of mod- 
els, and clearly deserves further study in the future. 
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